home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
CD Exchange
/
CD Exchange - Volume 1.iso
/
education
/
airfoil
/
airfoil.c
< prev
next >
Wrap
C/C++ Source or Header
|
1989-09-24
|
13KB
|
489 lines
/********************************************************************
* *
* Joukowski Airfoil Generator with Streamline and Pressure *
* Distribution Algorithms *
* *
* Written by: Russell Leighton 15 March 1987 *
* Lancaster, CA *
* Modified by: David Foster 19 June 1988 *
* Rochester, MI *
* To include effect of induced circulation *
* by a first order approximation. *
********************************************************************/
#include "airfoil.h"
main()
{
float rs,theta,h,vel,eta,S;
int psi;
ULONG MessageClass;
open_things();
do_about();
/****************/
/* Loop forever */
/****************/
for(;;)
{
while (Continue)
{
/****************************************************/
/* Wait, initially and after each plot, for user to */
/* bring up the double menu requester */
/****************************************************/
Wait(1L<<w->UserPort->mp_SigBit);
if((message = (struct IntuiMessage *)GetMsg(w->UserPort)) != NULL)
{
MessageClass = message->Class;
ReplyMsg(message);
if (MessageClass == REQVERIFY)
{
do_request();
break;
}
} /* end if */
} /* end while */
do_init();
if(mode)
{
/********************/
/* Plot Streamlines */
/********************/
FILL = TRUE;
for (psi = 12;psi > 0;--psi)
{
do_mess();
if(!Continue) break;
PENUP;
SetAPen(rp,(long)(psi+1));
SetBPen(rp,(long)(psi+1));
for (theta = 0.015;theta < TWO_PI;theta += PI/100)
{
vel = psi/(4.*velocity*r*sin(theta));
S = (1.+sin(alpha)/sin(theta));
vel = abs(S/(2.*vel));
eta = sqrt(vel*vel+1.0)-vel;
rs = r*(1.+eta)/(1.-eta);
transform(rs,theta);
} /* end for */
AreaEnd(rp);
} /* end for */
/*******************************/
/* Plot Stagnation Streamlines */
/*******************************/
do_mess();
if(Continue)
{
FILL = FALSE;
PENUP;
SetAPen(rp,1L);
h = (r-4.0)/40.0;
theta = -alpha;
for (rs = 4;rs >= r;rs += h)
{
transform(rs,theta);
}
PENUP;
theta = PI + alpha;
for (rs = 4;rs > r;rs += h)
{
transform(rs,theta);
}
} /* end if */
} /* end if */
else
{
/******************************/
/* Plot Pressure Distribution */
/******************************/
do_mess();
if(Continue)
{
FILL = TRUE;
PENUP;
SetAPen(rp,2L);
SetBPen(rp,2L);
for (theta = 0.0;theta <= TWO_PI;theta += PI/100)
{
rs = r+(sin(theta)+sin(alpha))*(sin(theta)+sin(alpha));
transform(rs,theta);
} /* end for */
AreaEnd(rp);
} /* end if */
} /* end else */
/****************/
/* Plot Airfoil */
/****************/
do_mess();
if(Continue)
{
FILL = TRUE;
PENUP;
rs = r;
SetAPen(rp,1L);
SetBPen(rp,1L);
for (theta = 0.0;theta <= TWO_PI;theta += PI/100)
{
transform(rs,theta);
}
AreaEnd(rp);
} /* end if */
} /* end forever */
} /* end main */
do_init()
{
float a0;
SetAPen(rp,0L);
RectFill(rp,(long)(XMIN+1),(long)(YMIN+1),(long)(XMAX-1),(long)(YMAX-1));
SetOPen(rp,1L);
/***********************************************************/
/* Calculate circle constants (circle center and radius) */
/* from airfoil constants through a reverse transformation */
/***********************************************************/
alpha = (float)angle*PI/180;
c = (float)camber/25;
t = (float)thickness/25;
a = 0;
b = c/2;
do
{
a0 = a;
a = t*(2*a0+1)/4/sqrt(b*b+2*a0+1);
b = c*(1+2*a)/(2+2*a);
}
while(abs(a-a0) > TOL);
r = sqrt(b*b+(a+1)*(a+1));
Continue = TRUE;
}
do_mess()
{
ULONG MessageClass;
/******************************************************************/
/* Check for double menu requester. This can be brought up at any */
/* time but can not be displayed unless the message is replied */
/* to, therefore this routine must be called periodically. */
/******************************************************************/
if((message = (struct IntuiMessage *)GetMsg(w->UserPort)) != NULL)
{
MessageClass = message->Class;
ReplyMsg(message);
if(MessageClass == REQVERIFY) do_request();
}
ixo = XCEN;
iyo = YCEN;
}
transform(rs,theta)
float rs,theta;
{
float x,y,z,u,v;
int PLOT;
long ix,iy,cx,cy;
/********************************************************************/
/* This is the Joukowski transformation routine. This is also where */
/* the plotting is done. Plotting is usually done using the area */
/* fill routines (AreaDraw & AreaMove). If FILL is true then the */
/* points are used to build up the area shape to be filled. See the */
/* Rom Kernal manual containing the graphics primatives for more */
/* details. If FILL is false then normal plotting is done. */
/********************************************************************/
x = rs*cos(theta-alpha)+a;
y = rs*sin(theta-alpha)+b;
z = 1/(x*x+y*y);
u = x*(1+z)*cos(alpha)-y*(1-z)*sin(alpha);
v = y*(1-z)*cos(alpha)+x*(1+z)*sin(alpha);
ix = (long)(SFAC*u+XCEN);
iy = (long)(YCEN-SFAC*v);
if(FILL)
{
PLOT = FALSE;
cx = ix;
cy = iy;
/*************************************************************/
/* This subroutine also clips the display area, hence the */
/* following statements. Contrary to the what the Rom Kernal */
/* manual states, all clipping must be done by the program */
/* if the area fill routines are used otherwise a nasty */
/* crash will result if plotting takes place out of bounds. */
/* Only the x values are clipped accurately since, for this */
/* routine only the x values at the boundary need to be */
/* accurate. The PEN parameter indicates the pen status for */
/* moves and draws as set by either PENUP or PENDOWN. */
/*************************************************************/
if((ix <= XMIN) && (ixo > XMIN))
{
cy += (float)(iyo-iy)*(XMIN-ix)/(ixo-ix);
cx = XMIN;
}
else if((ixo <= XMIN) && (ix > XMIN))
{
iyo += (float)(iy-iyo)*(XMIN-ixo)/(ix-ixo);
ixo = XMIN;
AreaDraw(rp,ixo,iyo);
PLOT = TRUE;
}
else if((ix >= XMAX) && (ixo < XMAX))
{
cy += (float)(iyo-iy)*(XMAX-ix)/(ixo-ix);
cx = XMAX;
}
else if((ixo >= XMAX) && (ix < XMAX))
{
iyo += (float)(iy-iyo)*(XMAX-ixo)/(ix-ixo);
ixo = XMAX;
AreaDraw(rp,ixo,iyo);
PLOT = TRUE;
}
if((ixo > XMIN) && (ixo < XMAX)) PLOT = TRUE;
if(cy < YMIN) cy = YMIN;
if(cy > YMAX) cy = YMAX;
if(PLOT)
{
if(PEN) AreaDraw(rp,cx,cy);
else { AreaMove(rp,cx,cy); PENDOWN; }
}
ixo = ix;
iyo = iy;
}
else
{
if(PEN) Draw(rp,ix,iy);
else { Move(rp,ix,iy); PENDOWN; }
}
}
do_request()
{